!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Semidiscretized ODE for the incompressible Navier-Stokes CG system
!!
!! \n
!!
!! This module represents the semidiscretized problem as an ODE, with
!! the layout required by the time integrators of \c
!! mod_time_integrators.
!!
!! Many of the code main variables and work arrays are defined as
!! module variables in the present module. This is done in order to
!! provide straightforward access to such variables to the ode
!! functions.
!<
module mod_cg_ns_ode

!-----------------------------------------------------------------------

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode

 use mod_sparse, only: &
   mod_sparse_initialized, &
   ! sparse types
   t_col,       &
   t_tri,       &
   t_pm_sk,     &
   ! overloaded operators
   matmul,      &
   transpose,   &
   clear

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_comm_world

 use mod_linsolver, only: &
   mod_linsolver_initialized, &
   c_linpb

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_testcases, only: &
   mod_testcases_initialized, &
   coeff_dir

 use mod_cg_ns_setup, only: &
   mod_cg_ns_setup_initialized, &
   ! miscellaneous
   max_char_len,   &
   ! Grid
   grid, bcs, &
   ddc_grid,  &
   ! FE space
   pressure_stabilization, add_div_lagmultiplier, &
   p_stab_type, p_stab_coeff,   &
   ubase,     pbase,     &
   udofs,     pdofs,     &
   ddc_udofs, ddc_pdofs, &
   ! Linear system
   t_nodal_field, t_linpb_mumps, linpb

 use mod_cg_ns_locmat, only: &
   mod_cg_ns_locmat_initialized, &
   t_bcs_error

 use mod_cg_ns_linsystem, only: &
   mod_cg_ns_linsystem_initialized, &
   t_ls,               &
   ns_mat, ns_packsol, &
   ns_unpacksol,       &
   ns_dirdata, ns_rhs

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cg_ns_ode_constructor, &
   mod_cg_ns_ode_destructor,  &
   mod_cg_ns_ode_initialized, &
   ! The following variables must be set by the main program
   implicit_advection, velocity_stabilization, u_stab_coeff, &
   antisym_adv, upwinding,                                   &
   mmmt, gn, dofs_nat, dofs_dir, gdofs_nat,                  &
   uuu1, uuu2, u_guess, fff, fff1, fff2, rhs1, linsys,       &
   a22ia21, a22ib2t, a22if2, a22i, a12a22i, b2a22i,          &
   ! Main ODE types and functions
   t_cg_ns_ode, new_cg_ns_ode, clear,                        &
   ! Diagnostics
   t_linsys, t_postproc

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 
 !> ODE type
 !!
 !! The state fields are allocatable components, to allow different
 !! variables to represent diffent states. The fields that represent
 !! the underlying discretization are module variables, so that they
 !! are implicitly shared by any solution state. An alternative
 !! approach would be using pointer fields in \c t_cg_ns_ode; we
 !! prefer here module variables to avoid a large number of pointer
 !! fields.
 type, extends(c_ode) :: t_cg_ns_ode
  real(wp), allocatable :: uuu(:)
 contains
  procedure, pass(tnd) :: rhs => cg_ns_ode_rhs ! not used
  procedure, pass(z)   :: add => cg_ns_ode_add
  procedure, pass(y)   :: stimes => cg_ns_ode_stimes
  procedure, pass(x)   :: solve  => cg_ns_ode_solve
 end type t_cg_ns_ode

 ! private members

 interface clear
   module procedure clear_cg_ns_ode
 end interface

! Module variables

 ! time stepping
 logical :: implicit_advection

 ! Local matrixes
 real(wp), allocatable :: a22ia21(:,:,:), a22ib2t(:,:,:), a22if2(:,:), &
   a22i(:,:,:), a12a22i(:,:,:), b2a22i(:,:,:)

 ! Advection
 logical :: antisym_adv

 ! Upwinding
 character(len=max_char_len) :: upwinding
 logical :: velocity_stabilization
 real(wp) :: u_stab_coeff

 ! Global matrix
 integer :: gn ! global system size, including domain decomposition
 real(wp), allocatable :: fff(:), fff1(:), fff2(:), uuu2(:), u_guess(:)
 real(wp), allocatable, target :: rhs1(:)
 type(t_nodal_field) :: uuu1
 type(t_pm_sk), save, target :: mmmt ! transposed matrix
 type(t_ls) :: linsys

 ! bcs variables
 integer, allocatable :: dofs_nat(:), dofs_dir(:)
 integer, allocatable, target :: gdofs_nat(:)

 ! timing variables: can be set by the ode subroutine to provide
 ! diagnostics to the main program
 real(t_realtime) :: t_linsys, t_postproc

 logical, protected ::               &
   mod_cg_ns_ode_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_cg_ns_ode'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_cg_ns_ode_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
      (mod_sparse_initialized.eqv..false.) .or. &
   (mod_mpi_utils_initialized.eqv..false.) .or. &
   (mod_linsolver_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
   (mod_testcases_initialized.eqv..false.) .or. &
 (mod_cg_ns_setup_initialized.eqv..false.) .or. &
(mod_cg_ns_locmat_initialized.eqv..false.) .or. &
(mod_cg_ns_linsystem_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cg_ns_ode_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_cg_ns_ode_initialized = .true.
 end subroutine mod_cg_ns_ode_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cg_ns_ode_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_cg_ns_ode_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_cg_ns_ode_initialized = .false.
 end subroutine mod_cg_ns_ode_destructor

!-----------------------------------------------------------------------
 
 pure subroutine new_cg_ns_ode(obj,n)
  integer, intent(in) :: n
  type(t_cg_ns_ode), intent(out) :: obj

   allocate( obj%uuu(n) )
   obj%uuu = -( 0.9_wp*huge(wp) ) ! set to undefined

 end subroutine new_cg_ns_ode
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_cg_ns_ode(obj)
  type(t_cg_ns_ode), intent(out) :: obj

   ! allocatable components are implicitly deallocated
 
 end subroutine clear_cg_ns_ode
 
!-----------------------------------------------------------------------

 subroutine cg_ns_ode_rhs(tnd,t,uuu)
  real(wp),            intent(in)    :: t   !< time level
  class(c_ode),        intent(in)    :: uuu !< present state
  class(t_cg_ns_ode), intent(inout) :: tnd !< tendency

   select type(uuu); type is(t_cg_ns_ode)

   tnd%uuu = -( 0.9_wp*huge(wp) )

   end select

 end subroutine cg_ns_ode_rhs

!-----------------------------------------------------------------------

 pure subroutine cg_ns_ode_add(z,x,y)
  class(c_ode),        intent(in)    :: x, y
  class(t_cg_ns_ode), intent(inout) :: z

   select type(x); type is(t_cg_ns_ode)
     select type(y); type is(t_cg_ns_ode)

   z%uuu = x%uuu + y%uuu

     end select
   end select

 end subroutine cg_ns_ode_add

!-----------------------------------------------------------------------

 pure subroutine cg_ns_ode_stimes(y,r,x)
  real(wp),            intent(in)    :: r
  class(c_ode),        intent(in)    :: x
  class(t_cg_ns_ode), intent(inout) :: y

   select type(x); type is(t_cg_ns_ode)

   y%uuu = r*x%uuu

   end select

 end subroutine cg_ns_ode_stimes

!-----------------------------------------------------------------------

 !> Implicit NS problem
 !!
 !! In this subroutine we recast the implicit problem of the nonlinear
 !! Navier-Stokes equations in the general form prescribed by the
 !! modules \c mod_time_integrators and \c mod_multistep. To this end,
 !! we use an affine rewriting
 !! \f{displaymath}{
 !!   f(t,u) = A(t,u)u + a(t,u)
 !! \f}
 !! of the continuous equation, where the precise form of \f$A\f$ and
 !! \f$a\f$ depends on whether the advective terms are treated
 !! explicitly or implicitly. The linear problem then reads (see also
 !! \ref constrained)
 !! \f{displaymath}{
 !!  \begin{array}{rcl}
 !!   \displaystyle u - \sigma A(t,u_l)u & = &
 !!   \displaystyle \left( b +\sigma a(t,u_l) \right)_u \\
 !!   \displaystyle B(t,u_l)u & = & 0
 !!  \end{array}
 !! \f}
 !! where \f$u\f$ is the complete vector of velocity and pressure
 !! degrees of freedom. Notice that the constraint \f$B\f$ can become
 !! nonlinear when using pressure stabilizations.
 !!
 !! The resulting system may allow the static condensation of certain
 !! velocity degrees of freedom: this is completely unrelated to the
 !! chosen time stepping. Details concerning the resulting linear
 !! system can be found in \c mod_cg_ns_linsystem.
 !<
 subroutine cg_ns_ode_solve(t,x,sigma,b,xl)
  real(wp),     intent(in) :: t, sigma
  class(c_ode), intent(in) :: b, xl
  class(t_cg_ns_ode), intent(inout) :: x

  type(t_bcs_error) :: b_err
  real(t_realtime) :: t0, t1
  character(len=*), parameter :: this_sub_name = 'cg_ns_ode_solve'

   select type(b); type is(t_cg_ns_ode)
     select type(xl); type is(t_cg_ns_ode)

   !--------------------------------------------------------------------
   ! Update system matrix and rhs
   !$ if(detailed_timing_omp) then
   !$   call omput_push_key("MatrixAndRhs")
   !$   call omput_start_timer()
   !$ endif
   if(.not.implicit_advection) then
   
     ! Right hand side
     call ns_rhs( fff,a22if2,linsys,ubase,grid,udofs,pdofs, &
                  b%uuu, xl%uuu, sigma, a22i,a12a22i,b2a22i )

     t_linsys = 0.0_t_realtime ! no time spent for the factorization
   else
   
     ! Build the matrix and the rhs
     !$ if(detailed_timing_omp) then
     !$   call omput_push_key("NsMat")
     !$   call omput_start_timer()
     !$ endif
     call ns_mat(mmmt,fff,linsys,a22ia21,a22ib2t,a22if2, ddc_grid , &
            ubase,pbase,grid,bcs,udofs,pdofs,add_div_lagmultiplier, &
   pressure_stabilization,trim(p_stab_type),velocity_stabilization, &
            b_err , bbb=b%uuu , uuu=xl%uuu , sigma=sigma )
     if(b_err%lerr) &
       call error(this_mod_name,this_sub_name,b_err%message)
     !$ if(detailed_timing_omp) then ! NsMat
     !$   call omput_write_time()
     !$   call omput_close_timer()
     !$   call omput_pop_key()
     !$ endif
     
     ! Factorize the system
     !$ if(detailed_timing_omp) then
     !$   call omput_push_key("factor")
     !$   call omput_start_timer()
     !$ endif
     t0 = my_second()
     call linpb%factor('factorization')
     t1 = my_second(); t_linsys = t1-t0
     !$ if(detailed_timing_omp) then ! factor
     !$   call omput_write_time()
     !$   call omput_close_timer()
     !$   call omput_pop_key()
     !$ endif

   endif
   !$ if(detailed_timing_omp) then ! MatrixAndRhs
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$ endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Implicit solve
   !$ if(detailed_timing_omp) then
   !$   call omput_push_key("ImplicitSolve")
   !$   call omput_start_timer()
   !$ endif
   t0 = my_second()

   !$ if(detailed_timing_omp) call omput_push_key("rhs")
   call ns_dirdata(uuu2,udofs,coeff_dir,ddc_udofs, &
                   velocity_stabilization)
   fff1 = fff(dofs_nat+1)
   rhs1 = fff1 - matmul(uuu2,mmmt%m(2,1))
   !$ if(detailed_timing_omp) call omput_pop_key()

   !$ if(detailed_timing_omp) call omput_push_key("solve")
   call ns_unpacksol(uuu1%x,xl%uuu,dofs_nat,linsys,grid)

   call linpb%solve(uuu1)
   !$ if(detailed_timing_omp) call omput_pop_key()

   ! boundary terms (diagnostics: compute fff2)
   !$ if(detailed_timing_omp) call omput_push_key("bcs")
   fff2 = matmul(uuu1%x,mmmt%m(1,2)) + matmul(uuu2,mmmt%m(2,2))
   fff(dofs_dir+1) = fff2
   !$ if(detailed_timing_omp) call omput_pop_key()

   t1 = my_second(); t_linsys = t_linsys + t1-t0
   !$ if(detailed_timing_omp) then ! ImplicitSolve
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$ endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Postprocessings
   !$ if(detailed_timing_omp) then
   !$   call omput_push_key("postprocessing")
   !$   call omput_start_timer()
   !$ endif
   t0 = my_second()

   call ns_packsol( x%uuu , ubase,pbase , grid,udofs,pdofs ,    &
     dofs_nat,dofs_dir,linsys,uuu1%x,uuu2,a22if2,a22ia21,a22ib2t)

   t1 = my_second(); t_postproc = t1-t0
   !$ if(detailed_timing_omp) then
   !$   call omput_write_time()
   !$   call omput_close_timer()
   !$   call omput_pop_key()
   !$ endif
   !--------------------------------------------------------------------

     end select
   end select

 end subroutine cg_ns_ode_solve

!-----------------------------------------------------------------------

end module mod_cg_ns_ode

